!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine tvd_mmc (ro, mx, my, mz, bx, by, bz, en, po, gx, gy, gz, x, y, z, dx, dy, dz, dt, nx, ny, nz,                           &
   fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd, fxpo_bnd,                                       &
   fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd, fypo_bnd,                                       &
   fzro_bnd, fzmx_bnd, fzmy_bnd, fzmz_bnd, fzbx_bnd, fzby_bnd, fzbz_bnd, fzen_bnd, fzpo_bnd)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny, nz

   real(8), intent(in) :: dx, dy, dz, dt, x(nx), y(ny), z(nz)
   real(8), dimension(nx, ny, nz), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, po, gx, gy, gz
   real(8), dimension(ny, nz, n_neighs), intent(inout) ::                                                                          &
      fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd, fxpo_bnd
   real(8), dimension(nx, nz, n_neighs), intent(inout) ::                                                                          &
      fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd, fypo_bnd
   real(8), dimension(ny, nz, n_neighs), intent(inout) ::                                                                          &
      fzro_bnd, fzmx_bnd, fzmy_bnd, fzmz_bnd, fzbx_bnd, fzby_bnd, fzbz_bnd, fzen_bnd, fzpo_bnd

   integer(4) :: ip1, im1, jp1, jm1, kp1, km1

   real(8) :: c2, s2, ca2, cfx, cfy, cfz
   real(8) :: ch, ch2, cp2
   real(8) :: dt_dx, dt_dy, dt_dz, dx2, dy2, dz2

   real(8), dimension(nx, ny, nz) :: fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen, fxpo
   real(8), dimension(nx, ny, nz) :: fyro, fymx, fymy, fymz, fybx, fyby, fybz, fyen, fypo
   real(8), dimension(nx, ny, nz) :: fzro, fzmx, fzmy, fzmz, fzbx, fzby, fzbz, fzen, fzpo
   real(8), dimension(nx, ny, nz) :: fxro_p, fxmx_p, fxmy_p, fxmz_p, fxbx_p, fxby_p, fxbz_p, fxen_p, fxpo_p
   real(8), dimension(nx, ny, nz) :: fxro_n, fxmx_n, fxmy_n, fxmz_n, fxbx_n, fxby_n, fxbz_n, fxen_n, fxpo_n
   real(8), dimension(nx, ny, nz) :: fyro_p, fymx_p, fymy_p, fymz_p, fybx_p, fyby_p, fybz_p, fyen_p, fypo_p
   real(8), dimension(nx, ny, nz) :: fyro_n, fymx_n, fymy_n, fymz_n, fybx_n, fyby_n, fybz_n, fyen_n, fypo_n
   real(8), dimension(nx, ny, nz) :: fzro_p, fzmx_p, fzmy_p, fzmz_p, fzbx_p, fzby_p, fzbz_p, fzen_p, fzpo_p
   real(8), dimension(nx, ny, nz) :: fzro_n, fzmx_n, fzmy_n, fzmz_n, fzbx_n, fzby_n, fzbz_n, fzen_n, fzpo_n
   real(8), dimension(nx, ny, nz) :: fxro_up, fxmx_up, fxmy_up, fxmz_up, fxbx_up, fxby_up, fxbz_up, fxen_up, fxpo_up
   real(8), dimension(nx, ny, nz) :: fyro_up, fymx_up, fymy_up, fymz_up, fybx_up, fyby_up, fybz_up, fyen_up, fypo_up
   real(8), dimension(nx, ny, nz) :: fzro_up, fzmx_up, fzmy_up, fzmz_up, fzbx_up, fzby_up, fzbz_up, fzen_up, fzpo_up
   real(8), dimension(nx, ny, nz) :: ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, po_pre

   real(8), dimension(nx, ny, nz) :: vx, vy, vz, v2, b2, pr, bx_resis, by_resis, bz_resis, en_resis, sx, sy, sz

   real(8), dimension(nx, ny, nz) :: fxro_pre_p, fxmx_pre_p, fxmy_pre_p, fxmz_pre_p,                                               &
      fxbx_pre_p, fxby_pre_p, fxbz_pre_p, fxen_pre_p, fxpo_pre_p
   real(8), dimension(nx, ny, nz) :: fxro_pre_n, fxmx_pre_n, fxmy_pre_n, fxmz_pre_n,                                               &
      fxbx_pre_n, fxby_pre_n, fxbz_pre_n, fxen_pre_n, fxpo_pre_n
   real(8), dimension(nx, ny, nz) :: fxro_lmt_p, fxmx_lmt_p, fxmy_lmt_p, fxmz_lmt_p,                                               &
      fxbx_lmt_p, fxby_lmt_p, fxbz_lmt_p, fxen_lmt_p, fxpo_lmt_p
   real(8), dimension(nx, ny, nz) :: fxro_lmt_n, fxmx_lmt_n, fxmy_lmt_n, fxmz_lmt_n,                                               &
      fxbx_lmt_n, fxby_lmt_n, fxbz_lmt_n, fxen_lmt_n, fxpo_lmt_n
   real(8), dimension(nx, ny, nz) :: fyro_pre_p, fymx_pre_p, fymy_pre_p, fymz_pre_p,                                               &
      fybx_pre_p, fyby_pre_p, fybz_pre_p, fyen_pre_p, fypo_pre_p
   real(8), dimension(nx, ny, nz) :: fyro_pre_n, fymx_pre_n, fymy_pre_n, fymz_pre_n,                                               &
      fybx_pre_n, fyby_pre_n, fybz_pre_n, fyen_pre_n, fypo_pre_n
   real(8), dimension(nx, ny, nz) :: fyro_lmt_p, fymx_lmt_p, fymy_lmt_p, fymz_lmt_p,                                               &
      fybx_lmt_p, fyby_lmt_p, fybz_lmt_p, fyen_lmt_p, fypo_lmt_p
   real(8), dimension(nx, ny, nz) :: fyro_lmt_n, fymx_lmt_n, fymy_lmt_n, fymz_lmt_n,                                               &
      fybx_lmt_n, fyby_lmt_n, fybz_lmt_n, fyen_lmt_n, fypo_lmt_n
   real(8), dimension(nx, ny, nz) :: fzro_pre_p, fzmx_pre_p, fzmy_pre_p, fzmz_pre_p,                                               &
      fzbx_pre_p, fzby_pre_p, fzbz_pre_p, fzen_pre_p, fzpo_pre_p
   real(8), dimension(nx, ny, nz) :: fzro_pre_n, fzmx_pre_n, fzmy_pre_n, fzmz_pre_n,                                               &
      fzbx_pre_n, fzby_pre_n, fzbz_pre_n, fzen_pre_n, fzpo_pre_n
   real(8), dimension(nx, ny, nz) :: fzro_lmt_p, fzmx_lmt_p, fzmy_lmt_p, fzmz_lmt_p,                                               &
      fzbx_lmt_p, fzby_lmt_p, fzbz_lmt_p, fzen_lmt_p, fzpo_lmt_p
   real(8), dimension(nx, ny, nz) :: fzro_lmt_n, fzmx_lmt_n, fzmy_lmt_n, fzmz_lmt_n,                                               &
      fzbx_lmt_n, fzby_lmt_n, fzbz_lmt_n, fzen_lmt_n, fzpo_lmt_n

!----------------------------------------------------------------------------------------------------------------------------------|
   dt_dx = dt / dx
   dt_dy = dt / dy
   dt_dz = dt / dz
   dx2 = dx * 2.0d0
   dy2 = dy * 2.0d0
   dz2 = dz * 2.0d0
   ch = cfl * minval ((/ dx, dy, dz /)) / min_dt * div_clean_flag
   ch2 = ch * ch
!   cp2 = ch * cd + eps
   cp2 = -dt * ch2 / log (cd) + eps

!----------------------------------------------------------------------------------------------------------------------------------|
!  step one
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx
      vx(i, j, k) = mx(i, j, k) / ro(i, j, k)
      vy(i, j, k) = my(i, j, k) / ro(i, j, k)
      vz(i, j, k) = mz(i, j, k) / ro(i, j, k)
      b2(i, j, k) = bx(i, j, k) * bx(i, j, k) + by(i, j, k) * by(i, j, k) + bz(i, j, k) * bz(i, j, k)
      v2(i, j, k) = vx(i, j, k) * vx(i, j, k) + vy(i, j, k) * vy(i, j, k) + vz(i, j, k) * vz(i, j, k)
      pr(i, j, k) = (en(i, j, k) - v2(i, j, k) * ro(i, j, k) / 2.0d0 - b2(i, j, k) / 2.0d0) * gmm1
      c2 = gm * pr(i, j, k)
      s2 = c2 + b2(i, j, k)
      ca2 = bx(i, j, k) * bx(i, j, k)
      cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro(i, j, k)  / 2.0d0)
      sx(i, j, k) = cfx + abs (vx(i, j, k))
      ca2 = by(i, j, k) * by(i, j, k)
      cfy = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro(i, j, k)  / 2.0d0)
      sy(i, j, k) = cfy + abs (vy(i, j, k))
      ca2 = bz(i, j, k) * bz(i, j, k)
      cfz = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro(i, j, k)  / 2.0d0)
      sz(i, j, k) = cfz + abs (vz(i, j, k))
      if (isnan(sx(i, j, k)) .or. sx(i, j, k) .gt. 1.0d0 / eps .or. sx(i, j, k) .lt. -1.0d0 / eps) sx(i, j, k) = 0.0d0
      if (isnan(sy(i, j, k)) .or. sy(i, j, k) .gt. 1.0d0 / eps .or. sy(i, j, k) .lt. -1.0d0 / eps) sy(i, j, k) = 0.0d0
      if (isnan(sz(i, j, k)) .or. sz(i, j, k) .gt. 1.0d0 / eps .or. sz(i, j, k) .lt. -1.0d0 / eps) sz(i, j, k) = 0.0d0
   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!     flux
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx
      fxro(i, j, k) = mx(i, j, k)
      fxmx(i, j, k) = vx(i, j, k) * mx(i, j, k) + pr(i, j, k) + b2(i, j, k) / 2.0d0 - bx(i, j, k) * bx(i, j, k)
      fxmy(i, j, k) = vx(i, j, k) * my(i, j, k) - bx(i, j, k) * by(i, j, k)
      fxmz(i, j, k) = vx(i, j, k) * mz(i, j, k) - bx(i, j, k) * bz(i, j, k)
      fxbx(i, j, k) = po(i, j, k)
      fxby(i, j, k) = vx(i, j, k) * by(i, j, k) - bx(i, j, k) * vy(i, j, k)
      fxbz(i, j, k) = vx(i, j, k) * bz(i, j, k) - bx(i, j, k) * vz(i, j, k)
      fxen(i, j, k) = (en(i, j, k) + pr(i, j, k) + b2(i, j, k) / 2.0d0) * vx(i, j, k) -                                            &
         (bx(i, j, k) * vx(i, j, k) + by(i, j, k) * vy(i, j, k) + bz(i, j, k) * vz(i, j, k)) * bx(i, j, k)
      fxpo(i, j, k) = ch2 * bx(i, j, k)

      fyro(i, j, k) = my(i, j, k)
      fymx(i, j, k) = vy(i, j, k) * mx(i, j, k) - by(i, j, k) * bx(i, j, k)
      fymy(i, j, k) = vy(i, j, k) * my(i, j, k) + pr(i, j, k) + b2(i, j, k) / 2.0d0 - by(i, j, k) * by(i, j, k)
      fymz(i, j, k) = vy(i, j, k) * mz(i, j, k) - by(i, j, k) * bz(i, j, k)
      fybx(i, j, k) = vy(i, j, k) * bx(i, j, k) - by(i, j, k) * vx(i, j, k)
      fyby(i, j, k) = po(i, j, k)
      fybz(i, j, k) = vy(i, j, k) * bz(i, j, k) - by(i, j, k) * vz(i, j, k)
      fyen(i, j, k) = (en(i, j, k) + pr(i, j, k) + b2(i, j, k) / 2.0d0) * vy(i, j, k) -                                            &
         (bx(i, j, k) * vx(i, j, k) + by(i, j, k) * vy(i, j, k) + bz(i, j, k) * vz(i, j, k)) * by(i, j, k)
      fypo(i, j, k) = ch2 * by(i, j, k)

      fzro(i, j, k) = mz(i, j, k)
      fzmx(i, j, k) = vz(i, j, k) * mx(i, j, k) - bz(i, j, k) * bx(i, j, k)
      fzmy(i, j, k) = vz(i, j, k) * my(i, j, k) - bz(i, j, k) * by(i, j, k)
      fzmz(i, j, k) = vz(i, j, k) * mz(i, j, k) + pr(i, j, k) + b2(i, j, k) / 2.0d0 - bz(i, j, k) * bz(i, j, k)
      fzbx(i, j, k) = vz(i, j, k) * bx(i, j, k) - bz(i, j, k) * vx(i, j, k)
      fzby(i, j, k) = vz(i, j, k) * by(i, j, k) - bz(i, j, k) * vy(i, j, k)
      fzbz(i, j, k) = po(i, j, k)
      fzen(i, j, k) = (en(i, j, k) + pr(i, j, k) + b2(i, j, k) / 2.0d0) * vz(i, j, k) -                                            &
         (bx(i, j, k) * vx(i, j, k) + by(i, j, k) * vy(i, j, k) + bz(i, j, k) * vz(i, j, k)) * bz(i, j, k)
      fzpo(i, j, k) = ch2 * bz(i, j, k)
      
!----------------------------------------------------------------------------------------------------------------------------------|
!     flux split
!----------------------------------------------------------------------------------------------------------------------------------|
      if (i .eq. 1) then
         max_speed = max (sx(i, j, k), sx(i + 1, j, k))
      elseif (i .eq. nx) then
         max_speed = max (sx(i - 1, j, k), sx(i, j, k))
      else
         max_speed = maxval (sx(i - 1 : i + 1, j, k))
      endif

      fxro_p(i, j, k) = (fxro(i, j, k) + max_speed * ro(i, j, k)) * 0.5d0
      fxro_n(i, j, k) = (fxro(i, j, k) - max_speed * ro(i, j, k)) * 0.5d0
      fxmx_p(i, j, k) = (fxmx(i, j, k) + max_speed * mx(i, j, k)) * 0.5d0
      fxmx_n(i, j, k) = (fxmx(i, j, k) - max_speed * mx(i, j, k)) * 0.5d0
      fxmy_p(i, j, k) = (fxmy(i, j, k) + max_speed * my(i, j, k)) * 0.5d0
      fxmy_n(i, j, k) = (fxmy(i, j, k) - max_speed * my(i, j, k)) * 0.5d0
      fxmz_p(i, j, k) = (fxmz(i, j, k) + max_speed * mz(i, j, k)) * 0.5d0
      fxmz_n(i, j, k) = (fxmz(i, j, k) - max_speed * mz(i, j, k)) * 0.5d0
      fxbx_p(i, j, k) = (fxbx(i, j, k) + max_speed * bx(i, j, k)) * 0.5d0
      fxbx_n(i, j, k) = (fxbx(i, j, k) - max_speed * bx(i, j, k)) * 0.5d0
      fxby_p(i, j, k) = (fxby(i, j, k) + max_speed * by(i, j, k)) * 0.5d0
      fxby_n(i, j, k) = (fxby(i, j, k) - max_speed * by(i, j, k)) * 0.5d0
      fxbz_p(i, j, k) = (fxbz(i, j, k) + max_speed * bz(i, j, k)) * 0.5d0
      fxbz_n(i, j, k) = (fxbz(i, j, k) - max_speed * bz(i, j, k)) * 0.5d0
      fxen_p(i, j, k) = (fxen(i, j, k) + max_speed * en(i, j, k)) * 0.5d0
      fxen_n(i, j, k) = (fxen(i, j, k) - max_speed * en(i, j, k)) * 0.5d0
      fxpo_p(i, j, k) = (fxpo(i, j, k) + max_speed * po(i, j, k)) * 0.5d0
      fxpo_n(i, j, k) = (fxpo(i, j, k) - max_speed * po(i, j, k)) * 0.5d0

      if (j .eq. 1) then
         max_speed = max (sy(i, j, k), sy(i, j + 1, k))
      elseif (j .eq. ny) then
         max_speed = max (sy(i, j - 1, k), sy(i, j, k))
      else
         max_speed = maxval (sy(i, j - 1 : j + 1, k))
      endif

      fyro_p(i, j, k) = (fyro(i, j, k) + max_speed * ro(i, j, k)) * 0.5d0
      fyro_n(i, j, k) = (fyro(i, j, k) - max_speed * ro(i, j, k)) * 0.5d0
      fymx_p(i, j, k) = (fymx(i, j, k) + max_speed * mx(i, j, k)) * 0.5d0
      fymx_n(i, j, k) = (fymx(i, j, k) - max_speed * mx(i, j, k)) * 0.5d0
      fymy_p(i, j, k) = (fymy(i, j, k) + max_speed * my(i, j, k)) * 0.5d0
      fymy_n(i, j, k) = (fymy(i, j, k) - max_speed * my(i, j, k)) * 0.5d0
      fymz_p(i, j, k) = (fymz(i, j, k) + max_speed * mz(i, j, k)) * 0.5d0
      fymz_n(i, j, k) = (fymz(i, j, k) - max_speed * mz(i, j, k)) * 0.5d0
      fybx_p(i, j, k) = (fybx(i, j, k) + max_speed * bx(i, j, k)) * 0.5d0
      fybx_n(i, j, k) = (fybx(i, j, k) - max_speed * bx(i, j, k)) * 0.5d0
      fyby_p(i, j, k) = (fyby(i, j, k) + max_speed * by(i, j, k)) * 0.5d0
      fyby_n(i, j, k) = (fyby(i, j, k) - max_speed * by(i, j, k)) * 0.5d0
      fybz_p(i, j, k) = (fybz(i, j, k) + max_speed * bz(i, j, k)) * 0.5d0
      fybz_n(i, j, k) = (fybz(i, j, k) - max_speed * bz(i, j, k)) * 0.5d0
      fyen_p(i, j, k) = (fyen(i, j, k) + max_speed * en(i, j, k)) * 0.5d0
      fyen_n(i, j, k) = (fyen(i, j, k) - max_speed * en(i, j, k)) * 0.5d0
      fypo_p(i, j, k) = (fypo(i, j, k) + max_speed * po(i, j, k)) * 0.5d0
      fypo_n(i, j, k) = (fypo(i, j, k) - max_speed * po(i, j, k)) * 0.5d0

      if (k .eq. 1) then
         max_speed = max (sz(i, j, k), sz(i, j, k + 1))
      elseif (k .eq. ny) then
         max_speed = max (sz(i, j, k - 1), sz(i, j, k))
      else
         max_speed = maxval (sz(i, j, k - 1 : k + 1))
      endif

      fzro_p(i, j, k) = (fzro(i, j, k) + max_speed * ro(i, j, k)) * 0.5d0
      fzro_n(i, j, k) = (fzro(i, j, k) - max_speed * ro(i, j, k)) * 0.5d0
      fzmx_p(i, j, k) = (fzmx(i, j, k) + max_speed * mx(i, j, k)) * 0.5d0
      fzmx_n(i, j, k) = (fzmx(i, j, k) - max_speed * mx(i, j, k)) * 0.5d0
      fzmy_p(i, j, k) = (fzmy(i, j, k) + max_speed * my(i, j, k)) * 0.5d0
      fzmy_n(i, j, k) = (fzmy(i, j, k) - max_speed * my(i, j, k)) * 0.5d0
      fzmz_p(i, j, k) = (fzmz(i, j, k) + max_speed * mz(i, j, k)) * 0.5d0
      fzmz_n(i, j, k) = (fzmz(i, j, k) - max_speed * mz(i, j, k)) * 0.5d0
      fzbx_p(i, j, k) = (fzbx(i, j, k) + max_speed * bx(i, j, k)) * 0.5d0
      fzbx_n(i, j, k) = (fzbx(i, j, k) - max_speed * bx(i, j, k)) * 0.5d0
      fzby_p(i, j, k) = (fzby(i, j, k) + max_speed * by(i, j, k)) * 0.5d0
      fzby_n(i, j, k) = (fzby(i, j, k) - max_speed * by(i, j, k)) * 0.5d0
      fzbz_p(i, j, k) = (fzbz(i, j, k) + max_speed * bz(i, j, k)) * 0.5d0
      fzbz_n(i, j, k) = (fzbz(i, j, k) - max_speed * bz(i, j, k)) * 0.5d0
      fzen_p(i, j, k) = (fzen(i, j, k) + max_speed * en(i, j, k)) * 0.5d0
      fzen_n(i, j, k) = (fzen(i, j, k) - max_speed * en(i, j, k)) * 0.5d0
      fzpo_p(i, j, k) = (fzpo(i, j, k) + max_speed * po(i, j, k)) * 0.5d0
      fzpo_n(i, j, k) = (fzpo(i, j, k) - max_speed * po(i, j, k)) * 0.5d0
   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  first-order upwind numerical flux
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz - 1
      kp1 = k + 1
      do j = 1, ny - 1
         jp1 = j + 1
         do i = 1, nx - 1
            ip1 = i + 1
            fxro_up(i, j, k) = fxro_p(i, j, k) + fxro_n(ip1, j, k)
            fxmx_up(i, j, k) = fxmx_p(i, j, k) + fxmx_n(ip1, j, k)
            fxmy_up(i, j, k) = fxmy_p(i, j, k) + fxmy_n(ip1, j, k)
            fxmz_up(i, j, k) = fxmz_p(i, j, k) + fxmz_n(ip1, j, k)
            fxbx_up(i, j, k) = fxbx_p(i, j, k) + fxbx_n(ip1, j, k)
            fxby_up(i, j, k) = fxby_p(i, j, k) + fxby_n(ip1, j, k)
            fxbz_up(i, j, k) = fxbz_p(i, j, k) + fxbz_n(ip1, j, k)
            fxen_up(i, j, k) = fxen_p(i, j, k) + fxen_n(ip1, j, k)
            fxpo_up(i, j, k) = fxpo_p(i, j, k) + fxpo_n(ip1, j, k)

            fyro_up(i, j, k) = fyro_p(i, j, k) + fyro_n(i, jp1, k)
            fymx_up(i, j, k) = fymx_p(i, j, k) + fymx_n(i, jp1, k)
            fymy_up(i, j, k) = fymy_p(i, j, k) + fymy_n(i, jp1, k)
            fymz_up(i, j, k) = fymz_p(i, j, k) + fymz_n(i, jp1, k)
            fybx_up(i, j, k) = fybx_p(i, j, k) + fybx_n(i, jp1, k)
            fyby_up(i, j, k) = fyby_p(i, j, k) + fyby_n(i, jp1, k)
            fybz_up(i, j, k) = fybz_p(i, j, k) + fybz_n(i, jp1, k)
            fyen_up(i, j, k) = fyen_p(i, j, k) + fyen_n(i, jp1, k)
            fypo_up(i, j, k) = fypo_p(i, j, k) + fypo_n(i, jp1, k)

            fzro_up(i, j, k) = fzro_p(i, j, k) + fzro_n(i, j, kp1)
            fzmx_up(i, j, k) = fzmx_p(i, j, k) + fzmx_n(i, j, kp1)
            fzmy_up(i, j, k) = fzmy_p(i, j, k) + fzmy_n(i, j, kp1)
            fzmz_up(i, j, k) = fzmz_p(i, j, k) + fzmz_n(i, j, kp1)
            fzbx_up(i, j, k) = fzbx_p(i, j, k) + fzbx_n(i, j, kp1)
            fzby_up(i, j, k) = fzby_p(i, j, k) + fzby_n(i, j, kp1)
            fzbz_up(i, j, k) = fzbz_p(i, j, k) + fzbz_n(i, j, kp1)
            fzen_up(i, j, k) = fzen_p(i, j, k) + fzen_n(i, j, kp1)
            fzpo_up(i, j, k) = fzpo_p(i, j, k) + fzpo_n(i, j, kp1)
         enddo
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  predictor-corrector step
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 2, nz
      km1 = k - 1
      do j = 2, ny
         jm1 = j - 1
         do i = 2, nx
            im1 = i - 1
            ro_pre(i, j, k) = ro(i, j, k) - dt_dx * (fxro_up(i, j, k) - fxro_up(im1, j, k)) -                                      &
                                            dt_dy * (fyro_up(i, j, k) - fyro_up(i, jm1, k)) -                                      &
                                            dt_dz * (fzro_up(i, j, k) - fzro_up(i, j, km1))
            mx_pre(i, j, k) = mx(i, j, k) - dt_dx * (fxmx_up(i, j, k) - fxmx_up(im1, j, k)) -                                      &
                                            dt_dy * (fymx_up(i, j, k) - fymx_up(i, jm1, k)) -                                      &
                                            dt_dz * (fzmx_up(i, j, k) - fzmx_up(i, j, km1))
            my_pre(i, j, k) = my(i, j, k) - dt_dx * (fxmy_up(i, j, k) - fxmy_up(im1, j, k)) -                                      &
                                            dt_dy * (fymy_up(i, j, k) - fymy_up(i, jm1, k)) -                                      &
                                            dt_dz * (fzmy_up(i, j, k) - fzmy_up(i, j, km1))
            mz_pre(i, j, k) = mz(i, j, k) - dt_dx * (fxmz_up(i, j, k) - fxmz_up(im1, j, k)) -                                      &
                                            dt_dy * (fymz_up(i, j, k) - fymz_up(i, jm1, k)) -                                      &
                                            dt_dz * (fzmz_up(i, j, k) - fzmz_up(i, j, km1))
            bx_pre(i, j, k) = bx(i, j, k) - dt_dx * (fxbx_up(i, j, k) - fxbx_up(im1, j, k)) -                                      &
                                            dt_dy * (fybx_up(i, j, k) - fybx_up(i, jm1, k)) -                                      &
                                            dt_dz * (fzbx_up(i, j, k) - fzbx_up(i, j, km1))
            by_pre(i, j, k) = by(i, j, k) - dt_dx * (fxby_up(i, j, k) - fxby_up(im1, j, k)) -                                      &
                                            dt_dy * (fyby_up(i, j, k) - fyby_up(i, jm1, k)) -                                      &
                                            dt_dz * (fzby_up(i, j, k) - fzby_up(i, j, km1))
            bz_pre(i, j, k) = bz(i, j, k) - dt_dx * (fxbz_up(i, j, k) - fxbz_up(im1, j, k)) -                                      &
                                            dt_dy * (fybz_up(i, j, k) - fybz_up(i, jm1, k)) -                                      &
                                            dt_dz * (fzbz_up(i, j, k) - fzbz_up(i, j, km1))
            en_pre(i, j, k) = en(i, j, k) - dt_dx * (fxen_up(i, j, k) - fxen_up(im1, j, k)) -                                      &
                                            dt_dy * (fyen_up(i, j, k) - fyen_up(i, jm1, k)) -                                      &
                                            dt_dz * (fzen_up(i, j, k) - fzen_up(i, j, km1))
            po_pre(i, j, k) = po(i, j, k) - dt_dx * (fxpo_up(i, j, k) - fxpo_up(im1, j, k)) -                                      &
                                            dt_dy * (fypo_up(i, j, k) - fypo_up(i, jm1, k)) -                                      &
                                            dt_dz * (fzpo_up(i, j, k) - fzpo_up(i, j, km1))
         enddo
      enddo
   enddo

   call resistivity(bx_resis, by_resis, bz_resis, en_resis, ro, mx, my, mz, bx, by, bz, en, x, y, z, dx, dy, dz, nx, ny, nz)

   call source_term(mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, po_pre,                                                &
      ro, mx, my, mz, bx, by, bz, po, gx, gy, gz, bx_resis, by_resis, bz_resis, en_resis,                                          &
      ch2, cp2, dx, dy, dz, dt, nx, ny, nz)

!----------------------------------------------------------------------------------------------------------------------------------|
!  step two
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx
!----------------------------------------------------------------------------------------------------------------------------------|
!     speed
!----------------------------------------------------------------------------------------------------------------------------------|
      vx(i, j, k) = mx_pre(i, j, k) / ro_pre(i, j, k)
      vy(i, j, k) = my_pre(i, j, k) / ro_pre(i, j, k)
      vz(i, j, k) = mz_pre(i, j, k) / ro_pre(i, j, k)
      b2(i, j, k) = bx_pre(i, j, k) * bx_pre(i, j, k) + by_pre(i, j, k) * by_pre(i, j, k) + bz_pre(i, j, k) * bz_pre(i, j, k)
      v2(i, j, k) = vx(i, j, k) * vx(i, j, k) + vy(i, j, k) * vy(i, j, k) + vz(i, j, k) * vz(i, j, k)
      pr(i, j, k) = (en_pre(i, j, k) - v2(i, j, k) * ro_pre(i, j, k) / 2.0d0 - b2(i, j, k) / 2.0d0) * gmm1
      c2 = gm * pr(i, j, k)
      s2 = c2 + b2(i, j, k)
      ca2 = bx_pre(i, j, k) * bx_pre(i, j, k)
      cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_pre(i, j, k) / 2.0d0)
      sx(i, j, k) = cfx + abs (vx(i, j, k))
      ca2 = by(i, j, k) * by(i, j, k)
      cfy = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_pre(i, j, k) / 2.0d0)
      sy(i, j, k) = cfy + abs (vy(i, j, k))
      ca2 = bz(i, j, k) * bz(i, j, k)
      cfz = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_pre(i, j, k) / 2.0d0)
      sz(i, j, k) = cfz + abs (vz(i, j, k))
      if (isnan(sx(i, j, k)) .or. sx(i, j, k) .gt. 1.0d0 / eps .or. sx(i, j, k) .lt. -1.0d0 / eps) sx(i, j, k) = 0.0d0
      if (isnan(sy(i, j, k)) .or. sy(i, j, k) .gt. 1.0d0 / eps .or. sy(i, j, k) .lt. -1.0d0 / eps) sy(i, j, k) = 0.0d0
      if (isnan(sz(i, j, k)) .or. sz(i, j, k) .gt. 1.0d0 / eps .or. sz(i, j, k) .lt. -1.0d0 / eps) sz(i, j, k) = 0.0d0
   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  flux
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx
      fxro(i, j, k) = mx_pre(i, j, k)
      fxmx(i, j, k) = vx(i, j, k) * mx_pre(i, j, k) + pr(i, j, k) + b2(i, j, k) / 2.0d0 - bx_pre(i, j, k) * bx_pre(i, j, k)
      fxmy(i, j, k) = vx(i, j, k) * my_pre(i, j, k) - bx_pre(i, j, k) * by_pre(i, j, k)
      fxmz(i, j, k) = vx(i, j, k) * mz_pre(i, j, k) - bx_pre(i, j, k) * bz_pre(i, j, k)
      fxbx(i, j, k) = po_pre(i, j, k)
      fxby(i, j, k) = vx(i, j, k) * by_pre(i, j, k) - bx_pre(i, j, k) * vy(i, j, k)
      fxbz(i, j, k) = vx(i, j, k) * bz_pre(i, j, k) - bx_pre(i, j, k) * vz(i, j, k)
      fxen(i, j, k) = (en_pre(i, j, k) + pr(i, j, k) + b2(i, j, k) / 2.0d0) * vx(i, j, k) -                                        &
         (bx_pre(i, j, k) * vx(i, j, k) + by_pre(i, j, k) * vy(i, j, k) + bz_pre(i, j, k) * vz(i, j, k)) * bx_pre(i, j, k)
      fxpo(i, j, k) = ch2 * bx_pre(i, j, k)

      fyro(i, j, k) = my_pre(i, j, k)
      fymx(i, j, k) = vy(i, j, k) * mx_pre(i, j, k) - by_pre(i, j, k) * bx_pre(i, j, k)
      fymy(i, j, k) = vy(i, j, k) * my_pre(i, j, k) + pr(i, j, k) + b2(i, j, k) / 2.0d0 - by_pre(i, j, k) * by_pre(i, j, k)
      fymz(i, j, k) = vy(i, j, k) * mz_pre(i, j, k) - by_pre(i, j, k) * bz_pre(i, j, k)
      fybx(i, j, k) = vy(i, j, k) * bx_pre(i, j, k) - by_pre(i, j, k) * vx(i, j, k)
      fyby(i, j, k) = po_pre(i, j, k)
      fybz(i, j, k) = vy(i, j, k) * bz_pre(i, j, k) - by_pre(i, j, k) * vz(i, j, k)
      fyen(i, j, k) = (en_pre(i, j, k) + pr(i, j, k) + b2(i, j, k) / 2.0d0) * vy(i, j, k) -                                        &
         (bx_pre(i, j, k) * vx(i, j, k) + by_pre(i, j, k) * vy(i, j, k) + bz_pre(i, j, k) * vz(i, j, k)) * by_pre(i, j, k)
      fypo(i, j, k) = ch2 * by_pre(i, j, k)

      fzro(i, j, k) = mz_pre(i, j, k)
      fzmx(i, j, k) = vz(i, j, k) * mx_pre(i, j, k) - bz_pre(i, j, k) * bx_pre(i, j, k)
      fzmy(i, j, k) = vz(i, j, k) * my_pre(i, j, k) - bz_pre(i, j, k) * by_pre(i, j, k)
      fzmz(i, j, k) = vz(i, j, k) * mz_pre(i, j, k) + pr(i, j, k) + b2(i, j, k) / 2.0d0 - bz_pre(i, j, k) * bz_pre(i, j, k)
      fzbx(i, j, k) = vz(i, j, k) * bx_pre(i, j, k) - bz_pre(i, j, k) * vx(i, j, k)
      fzby(i, j, k) = vz(i, j, k) * by_pre(i, j, k) - bz_pre(i, j, k) * vy(i, j, k)
      fzbz(i, j, k) = po_pre(i, j, k)
      fzen(i, j, k) = (en_pre(i, j, k) + pr(i, j, k) + b2(i, j, k) / 2.0d0) * vz(i, j, k) -                                        &
         (bx_pre(i, j, k) * vx(i, j, k) + by_pre(i, j, k) * vy(i, j, k) + bz_pre(i, j, k) * vz(i, j, k)) * bz_pre(i, j, k)
      fzpo(i, j, k) = ch2 * bz_pre(i, j, k)

!----------------------------------------------------------------------------------------------------------------------------------|
!     flux split
!----------------------------------------------------------------------------------------------------------------------------------|
      if (i .eq. 1) then
         max_speed = max (sx(i, j, k), sx(i + 1, j, k))
      elseif (i .eq. nx) then
         max_speed = max (sx(i - 1, j, k), sx(i, j, k))
      else
         max_speed = maxval (sx(i - 1 : i + 1, j, k))
      endif

      fxro_pre_p(i, j, k) = (fxro(i, j, k) + max_speed * ro_pre(i, j, k)) * 0.5d0
      fxro_pre_n(i, j, k) = (fxro(i, j, k) - max_speed * ro_pre(i, j, k)) * 0.5d0
      fxmx_pre_p(i, j, k) = (fxmx(i, j, k) + max_speed * mx_pre(i, j, k)) * 0.5d0
      fxmx_pre_n(i, j, k) = (fxmx(i, j, k) - max_speed * mx_pre(i, j, k)) * 0.5d0
      fxmy_pre_p(i, j, k) = (fxmy(i, j, k) + max_speed * my_pre(i, j, k)) * 0.5d0
      fxmy_pre_n(i, j, k) = (fxmy(i, j, k) - max_speed * my_pre(i, j, k)) * 0.5d0
      fxmz_pre_p(i, j, k) = (fxmz(i, j, k) + max_speed * mz_pre(i, j, k)) * 0.5d0
      fxmz_pre_n(i, j, k) = (fxmz(i, j, k) - max_speed * mz_pre(i, j, k)) * 0.5d0
      fxbx_pre_p(i, j, k) = (fxbx(i, j, k) + max_speed * bx_pre(i, j, k)) * 0.5d0
      fxbx_pre_n(i, j, k) = (fxbx(i, j, k) - max_speed * bx_pre(i, j, k)) * 0.5d0
      fxby_pre_p(i, j, k) = (fxby(i, j, k) + max_speed * by_pre(i, j, k)) * 0.5d0
      fxby_pre_n(i, j, k) = (fxby(i, j, k) - max_speed * by_pre(i, j, k)) * 0.5d0
      fxbz_pre_p(i, j, k) = (fxbz(i, j, k) + max_speed * bz_pre(i, j, k)) * 0.5d0
      fxbz_pre_n(i, j, k) = (fxbz(i, j, k) - max_speed * bz_pre(i, j, k)) * 0.5d0
      fxen_pre_p(i, j, k) = (fxen(i, j, k) + max_speed * en_pre(i, j, k)) * 0.5d0
      fxen_pre_n(i, j, k) = (fxen(i, j, k) - max_speed * en_pre(i, j, k)) * 0.5d0
      fxpo_pre_p(i, j, k) = (fxpo(i, j, k) + max_speed * po_pre(i, j, k)) * 0.5d0
      fxpo_pre_n(i, j, k) = (fxpo(i, j, k) - max_speed * po_pre(i, j, k)) * 0.5d0

      if (j .eq. 1) then
         max_speed = max (sy(i, j, k), sy(i, j + 1, k))
      elseif (j .eq. ny) then
         max_speed = max (sy(i, j - 1, k), sy(i, j, k))
      else
         max_speed = maxval (sy(i, j - 1 : j + 1, k))
      endif
      fyro_pre_p(i, j, k) = (fyro(i, j, k) + max_speed * ro_pre(i, j, k)) * 0.5d0
      fyro_pre_n(i, j, k) = (fyro(i, j, k) - max_speed * ro_pre(i, j, k)) * 0.5d0
      fymx_pre_p(i, j, k) = (fymx(i, j, k) + max_speed * mx_pre(i, j, k)) * 0.5d0
      fymx_pre_n(i, j, k) = (fymx(i, j, k) - max_speed * mx_pre(i, j, k)) * 0.5d0
      fymy_pre_p(i, j, k) = (fymy(i, j, k) + max_speed * my_pre(i, j, k)) * 0.5d0
      fymy_pre_n(i, j, k) = (fymy(i, j, k) - max_speed * my_pre(i, j, k)) * 0.5d0
      fymz_pre_p(i, j, k) = (fymz(i, j, k) + max_speed * mz_pre(i, j, k)) * 0.5d0
      fymz_pre_n(i, j, k) = (fymz(i, j, k) - max_speed * mz_pre(i, j, k)) * 0.5d0
      fybx_pre_p(i, j, k) = (fybx(i, j, k) + max_speed * bx_pre(i, j, k)) * 0.5d0
      fybx_pre_n(i, j, k) = (fybx(i, j, k) - max_speed * bx_pre(i, j, k)) * 0.5d0
      fyby_pre_p(i, j, k) = (fyby(i, j, k) + max_speed * by_pre(i, j, k)) * 0.5d0
      fyby_pre_n(i, j, k) = (fyby(i, j, k) - max_speed * by_pre(i, j, k)) * 0.5d0
      fybz_pre_p(i, j, k) = (fybz(i, j, k) + max_speed * bz_pre(i, j, k)) * 0.5d0
      fybz_pre_n(i, j, k) = (fybz(i, j, k) - max_speed * bz_pre(i, j, k)) * 0.5d0
      fyen_pre_p(i, j, k) = (fyen(i, j, k) + max_speed * en_pre(i, j, k)) * 0.5d0
      fyen_pre_n(i, j, k) = (fyen(i, j, k) - max_speed * en_pre(i, j, k)) * 0.5d0
      fypo_pre_p(i, j, k) = (fypo(i, j, k) + max_speed * po_pre(i, j, k)) * 0.5d0
      fypo_pre_n(i, j, k) = (fypo(i, j, k) - max_speed * po_pre(i, j, k)) * 0.5d0

      if (k .eq. 1) then
         max_speed = max (sz(i, j, k), sz(i, j, k + 1))
      elseif (k .eq. nz) then
         max_speed = max (sz(i, j, k - 1), sz(i, j, k))
      else
         max_speed = maxval (sz(i, j, k - 1 : k + 1))
      endif
      fzro_pre_p(i, j, k) = (fzro(i, j, k) + max_speed * ro_pre(i, j, k)) * 0.5d0
      fzro_pre_n(i, j, k) = (fzro(i, j, k) - max_speed * ro_pre(i, j, k)) * 0.5d0
      fzmx_pre_p(i, j, k) = (fzmx(i, j, k) + max_speed * mx_pre(i, j, k)) * 0.5d0
      fzmx_pre_n(i, j, k) = (fzmx(i, j, k) - max_speed * mx_pre(i, j, k)) * 0.5d0
      fzmy_pre_p(i, j, k) = (fzmy(i, j, k) + max_speed * my_pre(i, j, k)) * 0.5d0
      fzmy_pre_n(i, j, k) = (fzmy(i, j, k) - max_speed * my_pre(i, j, k)) * 0.5d0
      fzmz_pre_p(i, j, k) = (fzmz(i, j, k) + max_speed * mz_pre(i, j, k)) * 0.5d0
      fzmz_pre_n(i, j, k) = (fzmz(i, j, k) - max_speed * mz_pre(i, j, k)) * 0.5d0
      fzbx_pre_p(i, j, k) = (fzbx(i, j, k) + max_speed * bx_pre(i, j, k)) * 0.5d0
      fzbx_pre_n(i, j, k) = (fzbx(i, j, k) - max_speed * bx_pre(i, j, k)) * 0.5d0
      fzby_pre_p(i, j, k) = (fzby(i, j, k) + max_speed * by_pre(i, j, k)) * 0.5d0
      fzby_pre_n(i, j, k) = (fzby(i, j, k) - max_speed * by_pre(i, j, k)) * 0.5d0
      fzbz_pre_p(i, j, k) = (fzbz(i, j, k) + max_speed * bz_pre(i, j, k)) * 0.5d0
      fzbz_pre_n(i, j, k) = (fzbz(i, j, k) - max_speed * bz_pre(i, j, k)) * 0.5d0
      fzen_pre_p(i, j, k) = (fzen(i, j, k) + max_speed * en_pre(i, j, k)) * 0.5d0
      fzen_pre_n(i, j, k) = (fzen(i, j, k) - max_speed * en_pre(i, j, k)) * 0.5d0
      fzpo_pre_p(i, j, k) = (fzpo(i, j, k) + max_speed * po_pre(i, j, k)) * 0.5d0
      fzpo_pre_n(i, j, k) = (fzpo(i, j, k) - max_speed * po_pre(i, j, k)) * 0.5d0
   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  antidiffusive flux term
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz - 1
      kp1 = k + 1
      do j = 1, ny - 1
         jp1 = j + 1
         do i = 1, nx - 1
            ip1 = i + 1
            fxro_p(i, j, k) = fxro_pre_p(ip1, j, k) - fxro_p(i, j, k)
            fxro_n(i, j, k) = fxro_n(ip1, j, k) - fxro_pre_n(i, j, k)
            fxmx_p(i, j, k) = fxmx_pre_p(ip1, j, k) - fxmx_p(i, j, k)
            fxmx_n(i, j, k) = fxmx_n(ip1, j, k) - fxmx_pre_n(i, j, k)
            fxmy_p(i, j, k) = fxmy_pre_p(ip1, j, k) - fxmy_p(i, j, k)
            fxmy_n(i, j, k) = fxmy_n(ip1, j, k) - fxmy_pre_n(i, j, k)
            fxmz_p(i, j, k) = fxmz_pre_p(ip1, j, k) - fxmz_p(i, j, k)
            fxmz_n(i, j, k) = fxmz_n(ip1, j, k) - fxmz_pre_n(i, j, k)
            fxbx_p(i, j, k) = fxbx_pre_p(ip1, j, k) - fxbx_p(i, j, k)
            fxbx_n(i, j, k) = fxbx_n(ip1, j, k) - fxbx_pre_n(i, j, k)
            fxby_p(i, j, k) = fxby_pre_p(ip1, j, k) - fxby_p(i, j, k)
            fxby_n(i, j, k) = fxby_n(ip1, j, k) - fxby_pre_n(i, j, k)
            fxbz_p(i, j, k) = fxbz_pre_p(ip1, j, k) - fxbz_p(i, j, k)
            fxbz_n(i, j, k) = fxbz_n(ip1, j, k) - fxbz_pre_n(i, j, k)
            fxen_p(i, j, k) = fxen_pre_p(ip1, j, k) - fxen_p(i, j, k)
            fxen_n(i, j, k) = fxen_n(ip1, j, k) - fxen_pre_n(i, j, k)
            fxpo_p(i, j, k) = fxpo_pre_p(ip1, j, k) - fxpo_p(i, j, k)
            fxpo_n(i, j, k) = fxpo_n(ip1, j, k) - fxpo_pre_n(i, j, k)

            fyro_p(i, j, k) = fyro_pre_p(i, jp1, k) - fyro_p(i, j, k)
            fyro_n(i, j, k) = fyro_n(i, jp1, k) - fyro_pre_n(i, j, k)
            fymx_p(i, j, k) = fymx_pre_p(i, jp1, k) - fymx_p(i, j, k)
            fymx_n(i, j, k) = fymx_n(i, jp1, k) - fymx_pre_n(i, j, k)
            fymy_p(i, j, k) = fymy_pre_p(i, jp1, k) - fymy_p(i, j, k)
            fymy_n(i, j, k) = fymy_n(i, jp1, k) - fymy_pre_n(i, j, k)
            fymz_p(i, j, k) = fymz_pre_p(i, jp1, k) - fymz_p(i, j, k)
            fymz_n(i, j, k) = fymz_n(i, jp1, k) - fymz_pre_n(i, j, k)
            fybx_p(i, j, k) = fybx_pre_p(i, jp1, k) - fybx_p(i, j, k)
            fybx_n(i, j, k) = fybx_n(i, jp1, k) - fybx_pre_n(i, j, k)
            fyby_p(i, j, k) = fyby_pre_p(i, jp1, k) - fyby_p(i, j, k)
            fyby_n(i, j, k) = fyby_n(i, jp1, k) - fyby_pre_n(i, j, k)
            fybz_p(i, j, k) = fybz_pre_p(i, jp1, k) - fybz_p(i, j, k)
            fybz_n(i, j, k) = fybz_n(i, jp1, k) - fybz_pre_n(i, j, k)
            fyen_p(i, j, k) = fyen_pre_p(i, jp1, k) - fyen_p(i, j, k)
            fyen_n(i, j, k) = fyen_n(i, jp1, k) - fyen_pre_n(i, j, k)
            fypo_p(i, j, k) = fypo_pre_p(i, jp1, k) - fypo_p(i, j, k)
            fypo_n(i, j, k) = fypo_n(i, jp1, k) - fypo_pre_n(i, j, k)

            fzro_p(i, j, k) = fzro_pre_p(i, j, kp1) - fzro_p(i, j, k)
            fzro_n(i, j, k) = fzro_n(i, j, kp1) - fzro_pre_n(i, j, k)
            fzmx_p(i, j, k) = fzmx_pre_p(i, j, kp1) - fzmx_p(i, j, k)
            fzmx_n(i, j, k) = fzmx_n(i, j, kp1) - fzmx_pre_n(i, j, k)
            fzmy_p(i, j, k) = fzmy_pre_p(i, j, kp1) - fzmy_p(i, j, k)
            fzmy_n(i, j, k) = fzmy_n(i, j, kp1) - fzmy_pre_n(i, j, k)
            fzmz_p(i, j, k) = fzmz_pre_p(i, j, kp1) - fzmz_p(i, j, k)
            fzmz_n(i, j, k) = fzmz_n(i, j, kp1) - fzmz_pre_n(i, j, k)
            fzbx_p(i, j, k) = fzbx_pre_p(i, j, kp1) - fzbx_p(i, j, k)
            fzbx_n(i, j, k) = fzbx_n(i, j, kp1) - fzbx_pre_n(i, j, k)
            fzby_p(i, j, k) = fzby_pre_p(i, j, kp1) - fzby_p(i, j, k)
            fzby_n(i, j, k) = fzby_n(i, j, kp1) - fzby_pre_n(i, j, k)
            fzbz_p(i, j, k) = fzbz_pre_p(i, j, kp1) - fzbz_p(i, j, k)
            fzbz_n(i, j, k) = fzbz_n(i, j, kp1) - fzbz_pre_n(i, j, k)
            fzen_p(i, j, k) = fzen_pre_p(i, j, kp1) - fzen_p(i, j, k)
            fzen_n(i, j, k) = fzen_n(i, j, kp1) - fzen_pre_n(i, j, k)
            fzpo_p(i, j, k) = fzpo_pre_p(i, j, kp1) - fzpo_p(i, j, k)
            fzpo_n(i, j, k) = fzpo_n(i, j, kp1) - fzpo_pre_n(i, j, k)
         enddo
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  non-linear flux-limiter
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 2, nz
      km1 = k - 1
      do j = 2, ny
         jm1 = j - 1
         do i = 2, nx
            im1 = i - 1
            if (fxro_p(i, j, k) .eq. 0.0d0) fxro_p(i, j, k) = eps
            if (fxmx_p(i, j, k) .eq. 0.0d0) fxmx_p(i, j, k) = eps
            if (fxmy_p(i, j, k) .eq. 0.0d0) fxmy_p(i, j, k) = eps
            if (fxmz_p(i, j, k) .eq. 0.0d0) fxmz_p(i, j, k) = eps
            if (fxbx_p(i, j, k) .eq. 0.0d0) fxbx_p(i, j, k) = eps
            if (fxby_p(i, j, k) .eq. 0.0d0) fxby_p(i, j, k) = eps
            if (fxbz_p(i, j, k) .eq. 0.0d0) fxbz_p(i, j, k) = eps
            if (fxen_p(i, j, k) .eq. 0.0d0) fxen_p(i, j, k) = eps
            if (fxpo_p(i, j, k) .eq. 0.0d0) fxpo_p(i, j, k) = eps

            fxro_lmt_p(i, j, k) = fxro_p(im1, j, k) / fxro_p(i, j, k)
            fxmx_lmt_p(i, j, k) = fxmx_p(im1, j, k) / fxmx_p(i, j, k)
            fxmy_lmt_p(i, j, k) = fxmy_p(im1, j, k) / fxmy_p(i, j, k)
            fxmz_lmt_p(i, j, k) = fxmz_p(im1, j, k) / fxmz_p(i, j, k)
            fxbx_lmt_p(i, j, k) = fxbx_p(im1, j, k) / fxbx_p(i, j, k)
            fxby_lmt_p(i, j, k) = fxby_p(im1, j, k) / fxby_p(i, j, k)
            fxbz_lmt_p(i, j, k) = fxbz_p(im1, j, k) / fxbz_p(i, j, k)
            fxen_lmt_p(i, j, k) = fxen_p(im1, j, k) / fxen_p(i, j, k)
            fxpo_lmt_p(i, j, k) = fxpo_p(im1, j, k) / fxpo_p(i, j, k)

            fxro_lmt_p(i, j, k) = min(1.0d0, abs(fxro_lmt_p(i, j, k)))
            fxmx_lmt_p(i, j, k) = min(1.0d0, abs(fxmx_lmt_p(i, j, k)))
            fxmy_lmt_p(i, j, k) = min(1.0d0, abs(fxmy_lmt_p(i, j, k)))
            fxmz_lmt_p(i, j, k) = min(1.0d0, abs(fxmz_lmt_p(i, j, k)))
            fxbx_lmt_p(i, j, k) = min(1.0d0, abs(fxbx_lmt_p(i, j, k)))
            fxby_lmt_p(i, j, k) = min(1.0d0, abs(fxby_lmt_p(i, j, k)))
            fxbz_lmt_p(i, j, k) = min(1.0d0, abs(fxbz_lmt_p(i, j, k)))
            fxen_lmt_p(i, j, k) = min(1.0d0, abs(fxen_lmt_p(i, j, k)))
            fxpo_lmt_p(i, j, k) = min(1.0d0, abs(fxpo_lmt_p(i, j, k)))

            if (fyro_p(i, j, k) .eq. 0.0d0) fyro_p(i, j, k) = eps
            if (fymx_p(i, j, k) .eq. 0.0d0) fymx_p(i, j, k) = eps
            if (fymy_p(i, j, k) .eq. 0.0d0) fymy_p(i, j, k) = eps
            if (fymz_p(i, j, k) .eq. 0.0d0) fymz_p(i, j, k) = eps
            if (fybx_p(i, j, k) .eq. 0.0d0) fybx_p(i, j, k) = eps
            if (fyby_p(i, j, k) .eq. 0.0d0) fyby_p(i, j, k) = eps
            if (fybz_p(i, j, k) .eq. 0.0d0) fybz_p(i, j, k) = eps
            if (fyen_p(i, j, k) .eq. 0.0d0) fyen_p(i, j, k) = eps
            if (fypo_p(i, j, k) .eq. 0.0d0) fypo_p(i, j, k) = eps

            fyro_lmt_p(i, j, k) = fyro_p(i, jm1, k) / fyro_p(i, j, k)
            fymx_lmt_p(i, j, k) = fymx_p(i, jm1, k) / fymx_p(i, j, k)
            fymy_lmt_p(i, j, k) = fymy_p(i, jm1, k) / fymy_p(i, j, k)
            fymz_lmt_p(i, j, k) = fymz_p(i, jm1, k) / fymz_p(i, j, k)
            fybx_lmt_p(i, j, k) = fybx_p(i, jm1, k) / fybx_p(i, j, k)
            fyby_lmt_p(i, j, k) = fyby_p(i, jm1, k) / fyby_p(i, j, k)
            fybz_lmt_p(i, j, k) = fybz_p(i, jm1, k) / fybz_p(i, j, k)
            fyen_lmt_p(i, j, k) = fyen_p(i, jm1, k) / fyen_p(i, j, k)
            fypo_lmt_p(i, j, k) = fypo_p(i, jm1, k) / fypo_p(i, j, k)

            fyro_lmt_p(i, j, k) = min(1.0d0, abs(fyro_lmt_p(i, j, k)))
            fymx_lmt_p(i, j, k) = min(1.0d0, abs(fymx_lmt_p(i, j, k)))
            fymy_lmt_p(i, j, k) = min(1.0d0, abs(fymy_lmt_p(i, j, k)))
            fymz_lmt_p(i, j, k) = min(1.0d0, abs(fymz_lmt_p(i, j, k)))
            fybx_lmt_p(i, j, k) = min(1.0d0, abs(fybx_lmt_p(i, j, k)))
            fyby_lmt_p(i, j, k) = min(1.0d0, abs(fyby_lmt_p(i, j, k)))
            fybz_lmt_p(i, j, k) = min(1.0d0, abs(fybz_lmt_p(i, j, k)))
            fyen_lmt_p(i, j, k) = min(1.0d0, abs(fyen_lmt_p(i, j, k)))
            fypo_lmt_p(i, j, k) = min(1.0d0, abs(fypo_lmt_p(i, j, k)))

            if (fzro_p(i, j, k) .eq. 0.0d0) fzro_p(i, j, k) = eps
            if (fzmx_p(i, j, k) .eq. 0.0d0) fzmx_p(i, j, k) = eps
            if (fzmy_p(i, j, k) .eq. 0.0d0) fzmy_p(i, j, k) = eps
            if (fzmz_p(i, j, k) .eq. 0.0d0) fzmz_p(i, j, k) = eps
            if (fzbx_p(i, j, k) .eq. 0.0d0) fzbx_p(i, j, k) = eps
            if (fzby_p(i, j, k) .eq. 0.0d0) fzby_p(i, j, k) = eps
            if (fzbz_p(i, j, k) .eq. 0.0d0) fzbz_p(i, j, k) = eps
            if (fzen_p(i, j, k) .eq. 0.0d0) fzen_p(i, j, k) = eps
            if (fzpo_p(i, j, k) .eq. 0.0d0) fzpo_p(i, j, k) = eps

            fzro_lmt_p(i, j, k) = fzro_p(i, j, km1) / fzro_p(i, j, k)
            fzmx_lmt_p(i, j, k) = fzmx_p(i, j, km1) / fzmx_p(i, j, k)
            fzmy_lmt_p(i, j, k) = fzmy_p(i, j, km1) / fzmy_p(i, j, k)
            fzmz_lmt_p(i, j, k) = fzmz_p(i, j, km1) / fzmz_p(i, j, k)
            fzbx_lmt_p(i, j, k) = fzbx_p(i, j, km1) / fzbx_p(i, j, k)
            fzby_lmt_p(i, j, k) = fzby_p(i, j, km1) / fzby_p(i, j, k)
            fzbz_lmt_p(i, j, k) = fzbz_p(i, j, km1) / fzbz_p(i, j, k)
            fzen_lmt_p(i, j, k) = fzen_p(i, j, km1) / fzen_p(i, j, k)
            fzpo_lmt_p(i, j, k) = fzpo_p(i, j, km1) / fzpo_p(i, j, k)

            fzro_lmt_p(i, j, k) = min(1.0d0, abs(fzro_lmt_p(i, j, k)))
            fzmx_lmt_p(i, j, k) = min(1.0d0, abs(fzmx_lmt_p(i, j, k)))
            fzmy_lmt_p(i, j, k) = min(1.0d0, abs(fzmy_lmt_p(i, j, k)))
            fzmz_lmt_p(i, j, k) = min(1.0d0, abs(fzmz_lmt_p(i, j, k)))
            fzbx_lmt_p(i, j, k) = min(1.0d0, abs(fzbx_lmt_p(i, j, k)))
            fzby_lmt_p(i, j, k) = min(1.0d0, abs(fzby_lmt_p(i, j, k)))
            fzbz_lmt_p(i, j, k) = min(1.0d0, abs(fzbz_lmt_p(i, j, k)))
            fzen_lmt_p(i, j, k) = min(1.0d0, abs(fzen_lmt_p(i, j, k)))
            fzpo_lmt_p(i, j, k) = min(1.0d0, abs(fzpo_lmt_p(i, j, k)))
         enddo
      enddo
   enddo

   do k = 1, nz - 1
      kp1 = k + 1
      do j = 1, ny - 1
         jp1 = j + 1
         do i = 1, nx - 1
            ip1 = i + 1
            if (fxro_n(i, j, k) .eq. 0.0d0) fxro_n(i, j, k) = eps
            if (fxmx_n(i, j, k) .eq. 0.0d0) fxmx_n(i, j, k) = eps
            if (fxmy_n(i, j, k) .eq. 0.0d0) fxmy_n(i, j, k) = eps
            if (fxmz_n(i, j, k) .eq. 0.0d0) fxmz_n(i, j, k) = eps
            if (fxbx_n(i, j, k) .eq. 0.0d0) fxbx_n(i, j, k) = eps
            if (fxby_n(i, j, k) .eq. 0.0d0) fxby_n(i, j, k) = eps
            if (fxbz_n(i, j, k) .eq. 0.0d0) fxbz_n(i, j, k) = eps
            if (fxen_n(i, j, k) .eq. 0.0d0) fxen_n(i, j, k) = eps
            if (fxpo_n(i, j, k) .eq. 0.0d0) fxpo_n(i, j, k) = eps

            fxro_lmt_n(i, j, k) = fxro_n(ip1, j, k) / fxro_n(i, j, k)
            fxmx_lmt_n(i, j, k) = fxmx_n(ip1, j, k) / fxmx_n(i, j, k)
            fxmy_lmt_n(i, j, k) = fxmy_n(ip1, j, k) / fxmy_n(i, j, k)
            fxmz_lmt_n(i, j, k) = fxmz_n(ip1, j, k) / fxmz_n(i, j, k)
            fxbx_lmt_n(i, j, k) = fxbx_n(ip1, j, k) / fxbx_n(i, j, k)
            fxby_lmt_n(i, j, k) = fxby_n(ip1, j, k) / fxby_n(i, j, k)
            fxbz_lmt_n(i, j, k) = fxbz_n(ip1, j, k) / fxbz_n(i, j, k)
            fxen_lmt_n(i, j, k) = fxen_n(ip1, j, k) / fxen_n(i, j, k)
            fxpo_lmt_n(i, j, k) = fxpo_n(ip1, j, k) / fxpo_n(i, j, k)

            fxro_lmt_n(i, j, k) = min(1.0d0, abs(fxro_lmt_n(i, j, k)))
            fxmx_lmt_n(i, j, k) = min(1.0d0, abs(fxmx_lmt_n(i, j, k)))
            fxmy_lmt_n(i, j, k) = min(1.0d0, abs(fxmy_lmt_n(i, j, k)))
            fxmz_lmt_n(i, j, k) = min(1.0d0, abs(fxmz_lmt_n(i, j, k)))
            fxbx_lmt_n(i, j, k) = min(1.0d0, abs(fxbx_lmt_n(i, j, k)))
            fxby_lmt_n(i, j, k) = min(1.0d0, abs(fxby_lmt_n(i, j, k)))
            fxbz_lmt_n(i, j, k) = min(1.0d0, abs(fxbz_lmt_n(i, j, k)))
            fxen_lmt_n(i, j, k) = min(1.0d0, abs(fxen_lmt_n(i, j, k)))
            fxpo_lmt_n(i, j, k) = min(1.0d0, abs(fxpo_lmt_n(i, j, k)))

            if (fyro_n(i, j, k) .eq. 0.0d0) fyro_n(i, j, k) = eps
            if (fymx_n(i, j, k) .eq. 0.0d0) fymx_n(i, j, k) = eps
            if (fymy_n(i, j, k) .eq. 0.0d0) fymy_n(i, j, k) = eps
            if (fymz_n(i, j, k) .eq. 0.0d0) fymz_n(i, j, k) = eps
            if (fybx_n(i, j, k) .eq. 0.0d0) fybx_n(i, j, k) = eps
            if (fyby_n(i, j, k) .eq. 0.0d0) fyby_n(i, j, k) = eps
            if (fybz_n(i, j, k) .eq. 0.0d0) fybz_n(i, j, k) = eps
            if (fyen_n(i, j, k) .eq. 0.0d0) fyen_n(i, j, k) = eps
            if (fypo_n(i, j, k) .eq. 0.0d0) fypo_n(i, j, k) = eps

            fyro_lmt_n(i, j, k) = fyro_n(i, jp1, k) / fyro_n(i, j, k)
            fymx_lmt_n(i, j, k) = fymx_n(i, jp1, k) / fymx_n(i, j, k)
            fymy_lmt_n(i, j, k) = fymy_n(i, jp1, k) / fymy_n(i, j, k)
            fymz_lmt_n(i, j, k) = fymz_n(i, jp1, k) / fymz_n(i, j, k)
            fybx_lmt_n(i, j, k) = fybx_n(i, jp1, k) / fybx_n(i, j, k)
            fyby_lmt_n(i, j, k) = fyby_n(i, jp1, k) / fyby_n(i, j, k)
            fybz_lmt_n(i, j, k) = fybz_n(i, jp1, k) / fybz_n(i, j, k)
            fyen_lmt_n(i, j, k) = fyen_n(i, jp1, k) / fyen_n(i, j, k)
            fypo_lmt_n(i, j, k) = fypo_n(i, jp1, k) / fypo_n(i, j, k)

            fyro_lmt_n(i, j, k) = min(1.0d0, abs(fyro_lmt_n(i, j, k)))
            fymx_lmt_n(i, j, k) = min(1.0d0, abs(fymx_lmt_n(i, j, k)))
            fymy_lmt_n(i, j, k) = min(1.0d0, abs(fymy_lmt_n(i, j, k)))
            fymz_lmt_n(i, j, k) = min(1.0d0, abs(fymz_lmt_n(i, j, k)))
            fybx_lmt_n(i, j, k) = min(1.0d0, abs(fybx_lmt_n(i, j, k)))
            fyby_lmt_n(i, j, k) = min(1.0d0, abs(fyby_lmt_n(i, j, k)))
            fybz_lmt_n(i, j, k) = min(1.0d0, abs(fybz_lmt_n(i, j, k)))
            fyen_lmt_n(i, j, k) = min(1.0d0, abs(fyen_lmt_n(i, j, k)))
            fypo_lmt_n(i, j, k) = min(1.0d0, abs(fypo_lmt_n(i, j, k)))

            if (fzro_n(i, j, k) .eq. 0.0d0) fzro_n(i, j, k) = eps
            if (fzmx_n(i, j, k) .eq. 0.0d0) fzmx_n(i, j, k) = eps
            if (fzmy_n(i, j, k) .eq. 0.0d0) fzmy_n(i, j, k) = eps
            if (fzmz_n(i, j, k) .eq. 0.0d0) fzmz_n(i, j, k) = eps
            if (fzbx_n(i, j, k) .eq. 0.0d0) fzbx_n(i, j, k) = eps
            if (fzby_n(i, j, k) .eq. 0.0d0) fzby_n(i, j, k) = eps
            if (fzbz_n(i, j, k) .eq. 0.0d0) fzbz_n(i, j, k) = eps
            if (fzen_n(i, j, k) .eq. 0.0d0) fzen_n(i, j, k) = eps
            if (fzpo_n(i, j, k) .eq. 0.0d0) fzpo_n(i, j, k) = eps

            fzro_lmt_n(i, j, k) = fzro_n(i, j, kp1) / fzro_n(i, j, k)
            fzmx_lmt_n(i, j, k) = fzmx_n(i, j, kp1) / fzmx_n(i, j, k)
            fzmy_lmt_n(i, j, k) = fzmy_n(i, j, kp1) / fzmy_n(i, j, k)
            fzmz_lmt_n(i, j, k) = fzmz_n(i, j, kp1) / fzmz_n(i, j, k)
            fzbx_lmt_n(i, j, k) = fzbx_n(i, j, kp1) / fzbx_n(i, j, k)
            fzby_lmt_n(i, j, k) = fzby_n(i, j, kp1) / fzby_n(i, j, k)
            fzbz_lmt_n(i, j, k) = fzbz_n(i, j, kp1) / fzbz_n(i, j, k)
            fzen_lmt_n(i, j, k) = fzen_n(i, j, kp1) / fzen_n(i, j, k)
            fzpo_lmt_n(i, j, k) = fzpo_n(i, j, kp1) / fzpo_n(i, j, k)

            fzro_lmt_n(i, j, k) = min(1.0d0, abs(fzro_lmt_n(i, j, k)))
            fzmx_lmt_n(i, j, k) = min(1.0d0, abs(fzmx_lmt_n(i, j, k)))
            fzmy_lmt_n(i, j, k) = min(1.0d0, abs(fzmy_lmt_n(i, j, k)))
            fzmz_lmt_n(i, j, k) = min(1.0d0, abs(fzmz_lmt_n(i, j, k)))
            fzbx_lmt_n(i, j, k) = min(1.0d0, abs(fzbx_lmt_n(i, j, k)))
            fzby_lmt_n(i, j, k) = min(1.0d0, abs(fzby_lmt_n(i, j, k)))
            fzbz_lmt_n(i, j, k) = min(1.0d0, abs(fzbz_lmt_n(i, j, k)))
            fzen_lmt_n(i, j, k) = min(1.0d0, abs(fzen_lmt_n(i, j, k)))
            fzpo_lmt_n(i, j, k) = min(1.0d0, abs(fzpo_lmt_n(i, j, k)))
         enddo
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!   second-order numerical flux
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx
      fxro(i, j, k) = fxro_up(i, j, k) + (fxro_lmt_p(i, j, k) * fxro_p(i, j, k) - fxro_lmt_n(i, j, k) * fxro_n(i, j, k)) * 0.5d0
      fxmx(i, j, k) = fxmx_up(i, j, k) + (fxmx_lmt_p(i, j, k) * fxmx_p(i, j, k) - fxmx_lmt_n(i, j, k) * fxmx_n(i, j, k)) * 0.5d0
      fxmy(i, j, k) = fxmy_up(i, j, k) + (fxmy_lmt_p(i, j, k) * fxmy_p(i, j, k) - fxmy_lmt_n(i, j, k) * fxmy_n(i, j, k)) * 0.5d0
      fxmz(i, j, k) = fxmz_up(i, j, k) + (fxmz_lmt_p(i, j, k) * fxmz_p(i, j, k) - fxmz_lmt_n(i, j, k) * fxmz_n(i, j, k)) * 0.5d0
      fxbx(i, j, k) = fxbx_up(i, j, k) + (fxbx_lmt_p(i, j, k) * fxbx_p(i, j, k) - fxbx_lmt_n(i, j, k) * fxbx_n(i, j, k)) * 0.5d0
      fxby(i, j, k) = fxby_up(i, j, k) + (fxby_lmt_p(i, j, k) * fxby_p(i, j, k) - fxby_lmt_n(i, j, k) * fxby_n(i, j, k)) * 0.5d0
      fxbz(i, j, k) = fxbz_up(i, j, k) + (fxbz_lmt_p(i, j, k) * fxbz_p(i, j, k) - fxbz_lmt_n(i, j, k) * fxbz_n(i, j, k)) * 0.5d0
      fxen(i, j, k) = fxen_up(i, j, k) + (fxen_lmt_p(i, j, k) * fxen_p(i, j, k) - fxen_lmt_n(i, j, k) * fxen_n(i, j, k)) * 0.5d0
      fxpo(i, j, k) = fxpo_up(i, j, k) + (fxpo_lmt_p(i, j, k) * fxpo_p(i, j, k) - fxpo_lmt_n(i, j, k) * fxpo_n(i, j, k)) * 0.5d0

      fyro(i, j, k) = fyro_up(i, j, k) + (fyro_lmt_p(i, j, k) * fyro_p(i, j, k) - fyro_lmt_n(i, j, k) * fyro_n(i, j, k)) * 0.5d0
      fymx(i, j, k) = fymx_up(i, j, k) + (fymx_lmt_p(i, j, k) * fymx_p(i, j, k) - fymx_lmt_n(i, j, k) * fymx_n(i, j, k)) * 0.5d0
      fymy(i, j, k) = fymy_up(i, j, k) + (fymy_lmt_p(i, j, k) * fymy_p(i, j, k) - fymy_lmt_n(i, j, k) * fymy_n(i, j, k)) * 0.5d0
      fymz(i, j, k) = fymz_up(i, j, k) + (fymz_lmt_p(i, j, k) * fymz_p(i, j, k) - fymz_lmt_n(i, j, k) * fymz_n(i, j, k)) * 0.5d0
      fybx(i, j, k) = fybx_up(i, j, k) + (fybx_lmt_p(i, j, k) * fybx_p(i, j, k) - fybx_lmt_n(i, j, k) * fybx_n(i, j, k)) * 0.5d0
      fyby(i, j, k) = fyby_up(i, j, k) + (fyby_lmt_p(i, j, k) * fyby_p(i, j, k) - fyby_lmt_n(i, j, k) * fyby_n(i, j, k)) * 0.5d0
      fybz(i, j, k) = fybz_up(i, j, k) + (fybz_lmt_p(i, j, k) * fybz_p(i, j, k) - fybz_lmt_n(i, j, k) * fybz_n(i, j, k)) * 0.5d0
      fyen(i, j, k) = fyen_up(i, j, k) + (fyen_lmt_p(i, j, k) * fyen_p(i, j, k) - fyen_lmt_n(i, j, k) * fyen_n(i, j, k)) * 0.5d0
      fypo(i, j, k) = fypo_up(i, j, k) + (fypo_lmt_p(i, j, k) * fypo_p(i, j, k) - fypo_lmt_n(i, j, k) * fypo_n(i, j, k)) * 0.5d0

      fzro(i, j, k) = fzro_up(i, j, k) + (fzro_lmt_p(i, j, k) * fzro_p(i, j, k) - fzro_lmt_n(i, j, k) * fzro_n(i, j, k)) * 0.5d0
      fzmx(i, j, k) = fzmx_up(i, j, k) + (fzmx_lmt_p(i, j, k) * fzmx_p(i, j, k) - fzmx_lmt_n(i, j, k) * fzmx_n(i, j, k)) * 0.5d0
      fzmy(i, j, k) = fzmy_up(i, j, k) + (fzmy_lmt_p(i, j, k) * fzmy_p(i, j, k) - fzmy_lmt_n(i, j, k) * fzmy_n(i, j, k)) * 0.5d0
      fzmz(i, j, k) = fzmz_up(i, j, k) + (fzmz_lmt_p(i, j, k) * fzmz_p(i, j, k) - fzmz_lmt_n(i, j, k) * fzmz_n(i, j, k)) * 0.5d0
      fzbx(i, j, k) = fzbx_up(i, j, k) + (fzbx_lmt_p(i, j, k) * fzbx_p(i, j, k) - fzbx_lmt_n(i, j, k) * fzbx_n(i, j, k)) * 0.5d0
      fzby(i, j, k) = fzby_up(i, j, k) + (fzby_lmt_p(i, j, k) * fzby_p(i, j, k) - fzby_lmt_n(i, j, k) * fzby_n(i, j, k)) * 0.5d0
      fzbz(i, j, k) = fzbz_up(i, j, k) + (fzbz_lmt_p(i, j, k) * fzbz_p(i, j, k) - fzbz_lmt_n(i, j, k) * fzbz_n(i, j, k)) * 0.5d0
      fzen(i, j, k) = fzen_up(i, j, k) + (fzen_lmt_p(i, j, k) * fzen_p(i, j, k) - fzen_lmt_n(i, j, k) * fzen_n(i, j, k)) * 0.5d0
      fzpo(i, j, k) = fzpo_up(i, j, k) + (fzpo_lmt_p(i, j, k) * fzpo_p(i, j, k) - fzpo_lmt_n(i, j, k) * fzpo_n(i, j, k)) * 0.5d0
   enddo
   enddo
   enddo

   do k = 2, nz
      km1 = k - 1
      do j = 2, ny
         jm1 = j - 1
         do i = 2, nx
            im1 = i - 1
            ro(i, j, k) = ro(i, j, k) - dt_dx * (fxro(i, j, k) - fxro(im1, j, k)) -                                                &
                                        dt_dy * (fyro(i, j, k) - fyro(i, jm1, k)) -                                                &
                                        dt_dz * (fzro(i, j, k) - fzro(i, j, km1))
            mx(i, j, k) = mx(i, j, k) - dt_dx * (fxmx(i, j, k) - fxmx(im1, j, k)) -                                                &
                                        dt_dy * (fymx(i, j, k) - fymx(i, jm1, k)) -                                                &
                                        dt_dz * (fzmx(i, j, k) - fzmx(i, j, km1))
            my(i, j, k) = my(i, j, k) - dt_dx * (fxmy(i, j, k) - fxmy(im1, j, k)) -                                                &
                                        dt_dy * (fymy(i, j, k) - fymy(i, jm1, k)) -                                                &
                                        dt_dz * (fzmy(i, j, k) - fzmy(i, j, km1))
            mz(i, j, k) = mz(i, j, k) - dt_dx * (fxmz(i, j, k) - fxmz(im1, j, k)) -                                                &
                                        dt_dy * (fymz(i, j, k) - fymz(i, jm1, k)) -                                                &
                                        dt_dz * (fzmz(i, j, k) - fzmz(i, j, km1))
            bx(i, j, k) = bx(i, j, k) - dt_dx * (fxbx(i, j, k) - fxbx(im1, j, k)) -                                                &
                                        dt_dy * (fybx(i, j, k) - fybx(i, jm1, k)) -                                                &
                                        dt_dz * (fzbx(i, j, k) - fzbx(i, j, km1))
            by(i, j, k) = by(i, j, k) - dt_dx * (fxby(i, j, k) - fxby(im1, j, k)) -                                                &
                                        dt_dy * (fyby(i, j, k) - fyby(i, jm1, k)) -                                                &
                                        dt_dz * (fzby(i, j, k) - fzby(i, j, km1))
            bz(i, j, k) = bz(i, j, k) - dt_dx * (fxbz(i, j, k) - fxbz(im1, j, k)) -                                                &
                                        dt_dy * (fybz(i, j, k) - fybz(i, jm1, k)) -                                                &
                                        dt_dz * (fzbz(i, j, k) - fzbz(i, j, km1))
            en(i, j, k) = en(i, j, k) - dt_dx * (fxen(i, j, k) - fxen(im1, j, k)) -                                                &
                                        dt_dy * (fyen(i, j, k) - fyen(i, jm1, k)) -                                                &
                                        dt_dz * (fzen(i, j, k) - fzen(i, j, km1))
            po(i, j, k) = po(i, j, k) - dt_dx * (fxpo(i, j, k) - fxpo(im1, j, k)) -                                                &
                                        dt_dy * (fypo(i, j, k) - fypo(i, jm1, k)) -                                                &
                                        dt_dz * (fzpo(i, j, k) - fzpo(i, j, km1))
         enddo
      enddo
   enddo

   call resistivity(bx_resis, by_resis, bz_resis, en_resis,                                                                        &
      ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, x, y, z, dx, dy, dz, nx, ny, nz)

   call source_term(mx, my, mz, bx, by, bz, en, po,                                                                                &
      ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, po_pre, gx, gy, gz, bx_resis, by_resis, bz_resis, en_resis,          &
      ch2, cp2, dx, dy, dz, dt, nx, ny, nz)

!----------------------------------------------------------------------------------------------------------------------------------|
!  bnd value
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
      fxro_bnd(j, k, 1) = fxro(ng, j, k)
      fxro_bnd(j, k, 2) = fxro(nx - ng, j, k)

      fxmx_bnd(j, k, 1) = fxmx(ng, j, k)
      fxmx_bnd(j, k, 2) = fxmx(nx - ng, j, k)

      fxmy_bnd(j, k, 1) = fxmy(ng, j, k)
      fxmy_bnd(j, k, 2) = fxmy(nx - ng, j, k)

      fxmz_bnd(j, k, 1) = fxmz(ng, j, k)
      fxmz_bnd(j, k, 2) = fxmz(nx - ng, j, k)

      fxbx_bnd(j, k, 1) = fxbx(ng, j, k)
      fxbx_bnd(j, k, 2) = fxbx(nx - ng, j, k)

      fxby_bnd(j, k, 1) = fxby(ng, j, k)
      fxby_bnd(j, k, 2) = fxby(nx - ng, j, k)

      fxbz_bnd(j, k, 1) = fxbz(ng, j, k)
      fxbz_bnd(j, k, 2) = fxbz(nx - ng, j, k)

      fxen_bnd(j, k, 1) = fxen(ng, j, k)
      fxen_bnd(j, k, 2) = fxen(nx - ng, j, k)

      fxpo_bnd(j, k, 1) = fxpo(ng, j, k)
      fxpo_bnd(j, k, 2) = fxpo(nx - ng, j, k)
   enddo
   enddo

   do k = 1, nz
   do i = 1, nx
      fyro_bnd(i, k, 1) = fyro(i, ng, k)
      fyro_bnd(i, k, 2) = fyro(i, ny - ng, k)

      fymx_bnd(i, k, 1) = fymx(i, ng, k)
      fymx_bnd(i, k, 2) = fymx(i, ny - ng, k)

      fymy_bnd(i, k, 1) = fymy(i, ng, k)
      fymy_bnd(i, k, 2) = fymy(i, ny - ng, k)

      fymz_bnd(i, k, 1) = fymz(i, ng, k)
      fymz_bnd(i, k, 2) = fymz(i, ny - ng, k)

      fybx_bnd(i, k, 1) = fybx(i, ng, k)
      fybx_bnd(i, k, 2) = fybx(i, ny - ng, k)

      fyby_bnd(i, k, 1) = fyby(i, ng, k)
      fyby_bnd(i, k, 2) = fyby(i, ny - ng, k)

      fybz_bnd(i, k, 1) = fybz(i, ng, k)
      fybz_bnd(i, k, 2) = fybz(i, ny - ng, k)

      fyen_bnd(i, k, 1) = fyen(i, ng, k)
      fyen_bnd(i, k, 2) = fyen(i, ny - ng, k)

      fypo_bnd(i, k, 1) = fypo(i, ng, k)
      fypo_bnd(i, k, 2) = fypo(i, ny - ng, k)
   enddo
   enddo

   do j = 1, ny
   do i = 1, nx
      fzro_bnd(i, j, 1) = fzro(i, j, ng)
      fzro_bnd(i, j, 2) = fzro(i, j, nz - ng)

      fzmx_bnd(i, j, 1) = fzmx(i, j, ng)
      fzmx_bnd(i, j, 2) = fzmx(i, j, nz - ng)

      fzmy_bnd(i, j, 1) = fzmy(i, j, ng)
      fzmy_bnd(i, j, 2) = fzmy(i, j, nz - ng)

      fzmz_bnd(i, j, 1) = fzmz(i, j, ng)
      fzmz_bnd(i, j, 2) = fzmz(i, j, nz - ng)

      fzbx_bnd(i, j, 1) = fzbx(i, j, ng)
      fzbx_bnd(i, j, 2) = fzbx(i, j, nz - ng)

      fzby_bnd(i, j, 1) = fzby(i, j, ng)
      fzby_bnd(i, j, 2) = fzby(i, j, nz - ng)

      fzbz_bnd(i, j, 1) = fzbz(i, j, ng)
      fzbz_bnd(i, j, 2) = fzbz(i, j, nz - ng)

      fzen_bnd(i, j, 1) = fzen(i, j, ng)
      fzen_bnd(i, j, 2) = fzen(i, j, nz - ng)

      fzpo_bnd(i, j, 1) = fzpo(i, j, ng)
      fzpo_bnd(i, j, 2) = fzpo(i, j, nz - ng)
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine tvd_mmc
